With LinearSolve 4.0.0 supporting LinearProblem(A, B::Matrix) like A\B
(SciML/LinearSolve.jl#1072), PhiPadeCache now embeds a LinearSolve cache
and _phi_solve! runs the Pade denominator solve through it instead of
raw lu!/ldiv!. GenericLUFactorization is used because it refactorizes in
place with a cached pivot vector: workspace reuse is now fully
allocation-free (0 bytes measured for Float64 and ComplexF64 at
n = 7..100, p = 3..8, on Julia 1.10 and 1.12), improving on the lu!
version's per-call O(n) pivot allocation. (LUFactorization currently
copies A on every dense refactorization -- an in-place opt-in is being
added upstream, after which the LAPACK path can be swapped in for large
n.) Return-code semantics are unchanged: retcode Success -> info 0,
singular -> info 1 with NaN-filled outputs, non-finite results -> -1.

Tests tightened back to @allocated == 0; 254-assertion Phi testset
passes on 1.10 and 1.12 against LinearSolve 4.0.0.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* Allow LinearSolve v4

LinearSolve v4.0.0 adds batched (matrix) right-hand side support
(SciML/LinearSolve.jl#1072). NonlinearSolve is functionally unaffected —
the one production solve site in NonlinearSolveBase passes a flattened
vector — so this only widens compat to include v4 and bumps affected
subpackage patch versions. (NonlinearSolveBase already at 2.31.4.)

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

* Re-trigger CI against LinearSolve v4 (now registered in General)

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

* Allow LinearSolve v4

LinearSolve v4.0.0 adds batched (matrix) right-hand side support
(SciML/LinearSolve.jl#1072). OrdinaryDiffEq is functionally unaffected —
all call sites pass flattened vectors with explicit u0 — so this only
widens compat to include v4 and bumps affected subpackage patch versions.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

* Re-trigger CI against LinearSolve v4 (now registered in General)

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

…ix φ-functions (#236)

* Add O(p m^3) Al-Mohy--Liu scaling-and-recovering algorithm for phi functions

`phi(A, k)` / `phi!` previously computed phi_0(A), ..., phi_k(A) by calling
`phiv_dense!` on each of the m basis vectors, each exponentiating an
(m+k) x (m+k) block matrix -- overall O(m (m+k)^3).

This implements Algorithm 5.1 of Al-Mohy & Liu, "A scaling and recovering
algorithm for the matrix phi-functions" (arXiv:2506.01193, 2025), which
computes all phi functions simultaneously in O(k m^3): scale A by 2^-s,
evaluate the [m/m] diagonal Pade approximant to phi_p via Paterson--Stockmeyer,
recover the lower-index approximants with the recurrence
R^{(j)} = z R^{(j+1)} + 1/j!, and undo the scaling with the double-argument
formula. Backward-error-optimal Pade degree m and scaling s are selected from
the paper's theta table (Table 3.1) and the alpha_r / t cost analysis.

The new path is the default for Float64/ComplexF64 dense inputs (its Pade
tables are tuned for double precision); other element types and callers passing
the legacy `caches` bundle keep the basis-vector path. Diagonal inputs are
unchanged.

Validated against a high-precision block-exponential reference and the legacy
algorithm: worst-case relative error ~1e-14 across normal, non-normal,
large-norm, complex, Hessenberg and zero matrices for k = 1..10, matching or
beating the old accuracy. Measured speedups over the basis-vector path: 18x at
n=50, 51x at n=100, 72x at n=200 (k=3).

Closes #235.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* phi: reusable zero-alloc workspace + container-preserving static path

Address review feedback on #236:

* PhiPadeCache: a reusable workspace holding every scratch buffer. The
  strided Float64/ComplexF64 path (`_phi_almohy!`) now runs fully in place
  via mul!/getrf!/getrs! and broadcasting, allocating nothing once a cache
  exists (verified 0 bytes with @allocated across sizes). Pass it as the
  `caches` keyword to `phi!` to reuse across calls. (@stevengj: "do
  in-place, maybe?")

* Container-preserving generic path (`_phi_almohy_generic`) for immutable
  dense matrices: an SMatrix input stays an SMatrix throughout instead of
  being converted to a Matrix. Routed only for `!ismutable` inputs, so
  sparse/other mutable non-strided matrices keep the legacy path.
  (@stevengj: "unfortunate to use Matrix here if you start with SMatrix")

* Export and document PhiPadeCache and phi!; bump to 1.31.0 (new public
  API, SemVer minor).

Tests: workspace reuse + zero-allocation assertion, complex-matrix
workspace, and SMatrix static-in/static-out accuracy vs the block-exponential
reference. Full basictests + Aqua + JET pass locally (the only failures are
the pre-existing ForwardDiff-1.x Issue 42 asserts, green under CI's pinned
ForwardDiff 0.10.x).

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* phi: allocation-free LU across Julia versions (fixes LTS downgrade)

`LAPACK.getrf!(A, ipiv)` with a preallocated pivot vector only exists on
newer Julia; on the LTS (1.10) it errored, breaking the workspace solve.
ccall `getrf!` directly through libblastrampoline (as `StegrWork` does for
`stegr!`) so the factorization stays zero-allocation on all supported
versions; the 4-argument `getrs!` used for the solve is available
everywhere.

Verified on Julia 1.10 and 1.12: correctness ~2.5e-15 vs the BigFloat
reference and 0 bytes allocated on cache reuse.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* phi: ccall getrs! directly to keep the LU public-API-clean (QA)

The workspace solve used qualified `LinearAlgebra.BlasInt` and
`LinearAlgebra.LAPACK.getrs!`, which SciMLTesting's ExplicitImports
`all_qualified_accesses_are_public` check flags as reaching into non-public
internals (they are not on its ignore list). Mirror `exp_noalloc.jl`
instead: import `BlasInt` (an allow-listed explicit import) and ccall both
`getrf!` and `getrs!` through `BLAS.@blasfunc`/`BLAS.libblastrampoline`
(allow-listed names). No qualified access to non-public names remains in
`phi_almohy.jl`.

Behavior unchanged: verified correctness ~2.5e-15 and 0-byte cache reuse on
Julia 1.10 and 1.12; Aqua and JET still pass.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* phi: Horner PS, no-throw return codes, p-sized workspace, NaN guards

Re-analysis pass over the Al-Mohy--Liu implementation:

* Paterson--Stockmeyer rewritten in Horner form with the delta trick. The
  previous block scheme (mirroring the reference MATLAB) used explicit
  Atau powers plus two products per block -- up to 12 multiplications for
  m = 12 where Fasi's pi_m(tau) = 7, so the parameter selection was
  optimizing a cost the evaluator did not achieve. The Horner form's count
  now equals pi_m(tau) exactly. Measured 10-20% end-to-end speedup on the
  workspace path across n = 25..200.

* Numerical failure returns codes instead of throwing, so adaptive
  integrators can reject a step rather than abort: cache.info[] is 0 on
  success, the LAPACK getrf info for a singular Pade denominator, or -1
  for non-finite results (LAPACK does not flag NaN pivots, so a post-solve
  finiteness check backs the code); outputs are NaN-filled on failure.
  Fixes a latent InexactError where NaN/Inf input crashed in
  ceil(Int, t) during parameter selection before reaching the solve, and
  the generic path's SingularException from `\`. Wrongly sized workspaces
  (programmer errors) still throw, with informative messages.

* Workspace vectors sized from p at construction (was a fixed cap that
  BoundsError'd for p >~ 48) with an explicit capacity check; opnorm(A, 1)
  hoisted out of the per-degree ell loop (was recomputed 8x); dead Atau/Id
  buffers dropped (2n^2 less memory).

* Fix a closure-capture leak in the test helpers: phi_block_reference
  assigned to `n`, clobbering the enclosing testset's `n` (Julia closures
  assign to existing outer locals); the helpers now declare their
  temporaries `local`.

Verified: worst relative error 8.8e-15 over 33 BigFloat-referenced cases
(real/complex, normal/nonnormal, p = 1..8); 0 bytes allocated on cache
reuse (Float64 and ComplexF64, n = 10..100, p = 1..8, Julia 1.10 and
1.12); NaN/Inf inputs return info != 0 with NaN outputs and no throw;
full basictests 574 pass on 1.10 and 1.12 (only the pre-existing
ForwardDiff-1.x Issue 42 failures remain, which pass under CI's pinned
ForwardDiff 0.10); ExplicitImports clean.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* phi: overflow-safe backward-error coefficient; public lu! replaces ccalls

Two review-driven changes:

* _phi_be_coeff evaluated (m+p)!, m!, (2m+p)!, (2m+p+1)! separately before
  dividing, so intermediates could spuriously overflow (Inf/Inf = NaN for
  p >~ 159) even where the quotient is representable (@stevengj). Rewritten
  as a product of sub-unit factors,
  prod 1/(m+p+j) * prod 1/(m+j), which is monotonically decreasing and can
  only underflow when the true value does. Matches the BigFloat-exact value
  to 1.1e-15 over m = 1:12, p = 1:60; new p = 30 regression test.

* The hand-rolled getrf!/getrs! ccalls are gone: the solve now uses the
  public LinearAlgebra API, lu!(Dfact; check = false) + issuccess + ldiv!,
  which supports the matrix RHS directly and keeps the no-throw return-code
  semantics (issuccess false -> info > 0, NaN-filled outputs).
  LinearSolve.jl was considered but its LinearCache only solves vector
  right-hand sides through the public path, which would force a
  column-by-column loop or reaching into cache internals; plain lu!/ldiv!
  achieves ccall-free with zero new dependencies. Cost: the pivot vector is
  now the one per-call allocation (8n + ~100 bytes, measured 112 B at
  n = 7 to 928 B at n = 100; all O(n^2) buffers still reused), since no
  public API preallocates lu! pivots. Allocation tests assert that O(n)
  bound. The BlasInt import and ipiv field are removed; phi_almohy.jl is
  now entirely public-API.

Verified on Julia 1.10 and 1.12: 254-assertion Phi testset passes
(including new large-p and allocation-bound tests), full basictests 605
pass (only the pre-existing ForwardDiff-1.x Issue 42 failures remain),
ExplicitImports clean, NaN/Inf return codes intact.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* phi: solve through LinearSolve v4 batched matrix RHS

With LinearSolve 4.0.0 supporting LinearProblem(A, B::Matrix) like A\B
(SciML/LinearSolve.jl#1072), PhiPadeCache now embeds a LinearSolve cache
and _phi_solve! runs the Pade denominator solve through it instead of
raw lu!/ldiv!. GenericLUFactorization is used because it refactorizes in
place with a cached pivot vector: workspace reuse is now fully
allocation-free (0 bytes measured for Float64 and ComplexF64 at
n = 7..100, p = 3..8, on Julia 1.10 and 1.12), improving on the lu!
version's per-call O(n) pivot allocation. (LUFactorization currently
copies A on every dense refactorization -- an in-place opt-in is being
added upstream, after which the LAPACK path can be swapped in for large
n.) Return-code semantics are unchanged: retcode Success -> info 0,
singular -> info 1 with NaN-filled outputs, non-finite results -> -1.

Tests tightened back to @allocated == 0; 254-assertion Phi testset
passes on 1.10 and 1.12 against LinearSolve 4.0.0.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* Retrigger CI: LinearSolve 4.1.0 registered

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* exp_generic: fix dropped Dual partials at x = 0 under ForwardDiff >= 1

ForwardDiff 1.0 made iszero(::Dual) require zero partials as well, so
Dual(0, 1) fell through to intlog2, where ceil(Int, nx) dropped the
partials and 2^intlog2(0) = 2^64 overflowed to 0; the final
result^0 = one(result) then returned derivative 0 instead of 1 (Issue
42 testset). Gate the scaling on nx > 1 instead: values in (0, 1]
already produced s = 0, and Dual comparison uses the value only, so
zero-valued Duals with nonzero partials now take the s = 0 branch.
Verified against both ForwardDiff 0.10 and 1.x.

Widens the ForwardDiff test compat to "0.10.13, 1" accordingly; with
LinearSolve 4 in the dependency tree, PureKLU 1.1 forces ForwardDiff
>= 1, which is how this latent failure surfaced.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

---------

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>